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INTRODUCTION 


In  "Crack  Tip  Finite  Elements  are  Unnecessary",  Henshell  and 
Shaw  [1]  reported  that  the  inverse  of  the  Jacobian  associated  with  the 
coordinate  transformation  becomes  singular  at  a point  when  the  mid-side 
nodes  for  two-dimensional  eight-node  quadrilateral  elements  are  placed 
at  quarter  points.  Interestingly  enough,  the  same  singularity  was 
discovered  independently  by  Barsoum  [2]  for  two-dimensional,  as  well 
as  three-dimensional  quadratic  isoparametric  elements.  It  was  then 
natural  to  investigate  the  order  of  the  singularity,  and  it  was  found 
that  the  singularity  was  precisely  of  the  order  one-half  for  the 
strains,  a phenomenon  encountered  in  linear  fracture  mechanics.  This 
remarkable  phenomenon  completely  eliminates  the  necessity  of  incor- 
porating special  crack-tip  elements  [3,4,5]  and  has  additional  advan- 
tages over  the  special  crack-tip  elements;  namely,  it  satisfies 
constant  strain  and  rigid  body  modes.  The  special  crack-tip  elements 
were  introduced  in  the  literature  to  avoid  the  extremely  fine  grid 
mesh  required  in  the  vicinity  of  the  crack  and  the  cumbersome  extra- 
polation needed  when  using  regular  finite  elements  [6,7]. 

Advanced  versions  of  NASTRAN  [8],  as  well  as  some  general  purpose 
programs  have  such  isoparametric  elements.  Hence,  by  judicious  choice 
of  nodes,  accurate  crack-tip  elements  can  be  formulated  and  stress 
intensity  factors  for  cracks  and  flaws  can  be  computed. 

In  the  present  report,  after  a brief  review  of  the  theoretical 

[formulation,  we  discuss  the  implementation  of  these  elements  as  NASTRAN 
dummy  user  elements  [18].  Test  problems  are  done  to  assess  the 
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accuracy.  Stress  intensity  factors  are  computed  for  a C-shaped 
specimen. 


Finally  the  stability  of  two-dimensional  quadrilateral  element  as 
well  as  triangular  element  is  investigated;  and  it  is  shown  that  a 
triangular  element  preserves  the  required  singularity  under  perturbation 
of  the  mid-side  node.  [19] 

The  method  can  be  extended  to  higher  order  elements. 


SYMBOLS 

(x,y) 
U.  n) 


u,v 

{el 

[J] 

{a} 

CD] 

[K] 

E,G,v 

L L 3 

KrKn 

r,e 
{ F }e 


Cartesian  coordinates 
Curvilinear  coordinates 
Grid  point  coordinates 

Shape  function  at  grid  point  i 

Cartesian  displacements 
Strain  vector 
Jacobian  matrix 
Stress  vector 
Stress-strain  matrix 
Element  stiffness  matrix 
Elastic  constants 
Area  coordinates 
Stress  intensity  factors 
Local  cylindrical  coordinates 
Equivalent  nodal  forces 


THE  TWO-DIMENSIONAL  CASE 


QUADRILATERAL 

Following  the  notation  of  Zienkiewicz  [9]  the  eight  node  quadri- 
lateral element  in  Cartesian  coordinates  (x,y)  is  formulated  by  mapping 
its  geometry  into  the  curvilinear  space  (?,n)  of  the  normalized  square 
(-1  £ C 1 1 » -1  ± n ± 1 ) by  quadratic  shape  functions  of  the 
'Serendipity'  family  [9]: 

8 

x = .£  N. (c.n)x-  , 

1=  l 1 1 

8 

y =.t  N.(f,n)y.  , 

i=ii  i 

= [(l+CCi  )0+nn} ) - (1  -€2 ) (1+nn^ ) - (l-n2)(lH£i  )]t?n?/4 

(1) 

+ (l-C2)(l+nn« )(l-C?)n?/2  + (l-n2)(l+f f. )(1 -n2  K?/2, 

'll  'll 

where  N.  is  the  shape  function  at  node  i whose  Cartesian  and  curvilinear 
coordinates  are  (x^,y^)  and  ( C ^ , n ^ ) respectively,  (i  = l ...  8).  The 
details  of  the  shape  functions  and  the  numbering  sequence  are  given  in 
Figure  1.  The  same  shape  functions  are  used  to  interpolate  the  dis- 
placements within  the  element,  hence  the  name  isoparametric. 

0 

u =.£  N- (f,n)u-  , 

i=  l 1 1 

8 (2) 


■ 


3 


BMH 


lA(C-l)(n+l)U-n+l)  Ko  = i/U(i+c)(l*n)U+n-l) 


Figure  1.  Shape  functions  and  numbering  sequence  for  a quadrilateral  element. 


The  stiffness  matrix  is  found  in  the  usual  way  as  follows: 


e 

“a 

X 

3x  0 

u 

a G 

' y 

K 

0 3 

1 

» 

3y 

Y 

3 3_ 

V 

xy 

3y  3x 

• v, 

(3) 


Substituting  from  Equation  (2)  into  Equation  (3)  we  have: 

. *\ 

U.‘  I [ U; 

< 


{E}  = [B] 


V - [...  B.  ...] 


v. 


where 


[B.]  ■ 


?x 

0 

3Ni 

«y~ 

3Ni 

3Ni 

sy 

3X~ 

By  the  rules  of  partial  differentiation  we  obtain 


5Ni 

3Ni  ' 

3X 

■< 

- = [J]'1  - 

aT 

3Ni 

aNj 

W 

V.  -/ 

. - 

(a) 


(5) 


(6) 
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where  [J],  the  Jacobian  matrix,  by  virtue  of  Equation  (1)  is  given  by 

3N. 


[J]  i 


3X 

' 

35 

35 

3X 

3n 

3n 

35 


3^ 

3n 


xi  *i 


(7) 


The  stress  components  are  given  by 


(a)  =<a 


y=  [d]  { c j 


xy 


(8) 


where  [D]  is  the  stress-strain  matrix  and  for  the  case  of  plane  stress 
is  given  by 


1 v 0 


v 1 0 

0 0 (1 -v)/2 


The  element  stiffness  matrix  is  then: 


(9) 


1 1 

[K]  = f j*  [B]T[D][B]det  | J |d?,dn  . (10) 

-1  -1 

The  integration  in  Equation  (10)  is  done  numerically  by  nine-point 
Gaussian  quadrature  as  explained  in  [9]. 


..  ...  - . 
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TRIANGULAR 


Again,  following  the  notation  of  Zienkiewicz  [9]  the  six  node 
triangular  element  in  cartesian  coordinates  is  formulated  by  trans 
formations  using  area  coordinates  (L1.L2.L3): 


(11) 


(12) 

Equation  (12)  gives  the  shape  function  at  node  i whose  cartesian  and  area 
coordinates  are  (x^  ,y..)  and  (L^L^L^  respectively.  The  details  of  the 

shape  function  (12)  and  the  numbering  sequence  are  given  in  Figure  2. 

The  displacements  for  the  isoparametric  triangular  element  are  given  by: 

u = l N. (Lj ,L2,L3)ui  , 

V (13) 

v = £ (Lj ,L2,L3)v^  . 

1=  1 

The  rest  of  the  analysis  proceeds  in  a similar  way  as  for  the  quadrilat- 
eral case  given  above.  The  numerical  integration  for  the  stiffness 
matrix  is  done  using  Hammer's  Quadrature  [9]. 

7 


6 

x = £ N.j  (L  i,L2,L  3)x.j  , 
i=  1 

s 

y = £ N- (L  i,L2.L  3>y  , 

i=  1 1 

L 1 + L2  + L 3 = 1 , 

where 

",  * LlLu(2Ll-l)(2L11-l)  - 16L1L2LljL2l(2L3rl) 

* L2Laj(2L2-l)(2L2rl)  - 16L2L3L21L3l(2L1(-n 

* ‘-s<-31C2L3-1>(2L31-1)  - 16L3L,Ls1llj(2L2(-l)  , 
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THE  CRACK  TIP  SINGULARITY  ELEMENT 

It  is  clear  from  Equations  (4)  and  (6)  that  we  need  the  inverse  of 
Jacobian  matrix  [J]  before  the  strains  can  be  computed.  Hence,  whenever 
the  inverse  of  [J]  is  singular  or,  equivalently,  the  determinant  of  [J] 
is  zero,  the  strains  and  stresses  may  become  singular.  A singular 
Jacobian  is  achieved  by  placing  the  adjacent  mid-side  nodes  at  quarter 
points  from  the  node  where  the  singularity  is  desired. 

This  will  be  illustrated  by  investigating  the  singularity  at  node 
1 along  the  line  1 - 2 of  Figures  1 and  2. 

QUADRILATERAL  CASE 

Evaluating  the  shape  functions  given  in  Figure  1 at  n = -1 , we 
have  the  transformation 

x = -imO-OXj  + 1/2L(1  + C)x2  + (1-C2)x5  . (14) 

Choosing  Xj  = 0,  x2  = L and  the  quarter  point  x5  = L/4,  Fquation  (14) 
becomes 

x = 1/2CO+OL  + (l-C?)L/4  (15) 

Solving  for  £ we  have 

£ = -1  + 2/i7i  . (16) 

In  this  case  the  reduced  Jacobian  becomes 

||=  t(lH)  = */xU  . (17) 

Equation  (17)  clearly  indicates  the  singularity  for  the  inverse  of  the 
Jacobian  at  x = 0,  £ = -1.  The  order  of  singularity  can  be  obtained 
from  the  displacement  along  line  1 - 2 (Fig.  1).  From  Equation  (2) 
we  have 

u(£,-l ) = -l/2c(l-£)ui  + l/2f(l+£)u2  + (1-£2)u5  , 
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(18) 


and  writing  in  terms  of  x from  Equation  (16)  we  have 

u * -l/2(-H2^)(2-2^)uj  + 1/2(-1+2^)(2^)u2 

* 4(£-  ^u5  . 

Differentiating  Equation  (18)  we  obtain  the  strains  in  the  x-direction: 

'x  ■ If  • 'VM*?  - &«.  * + * C/xT  - 

indicating  the  singularity  of  order  one-half  (^),  precisely  the 

singularity  needed  for  crack  problems.  It  can  be  seen  that  Equation 
(19)  also  incorporates  constant  strain  terms. 

TRIANGULAR  CASE 

Evaluating  the  shape  functions  given  in  Figure  2 at  L3  = 0,  we 

have 

x = Lj^Lj-Dxj  + (1-L1)(1-2L1  )Xj  + 4L1(1-L1)x3  (20) 

choosing  x3  = 0,  x2  = A and  the  quarter  point  x4  = A/4,  Equation  (20) 


becomes 

x/A  = Lj  - 2Ll  + 1 (21) 

Solving  for  Li  we  have 

Li  = 1 - /x/A  (22) 

For  this  case  the  reduced  Jacobian  becomes 

|f  = 2A(Lj  -1 ) = -2^x  (23) 

3Li 

which  shows  the  singularity  for  the  inverse  of  the  Jacobian  at  x = 0, 
Lj  = 1.  The  order  of  singularity  can  be  obtained,  as  before,  from 
Equation  (13). 
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u(L3=0)  = (1-/x7A)(1-2^7A)Uj  + v/57A(2yx7A-l  )u2 


(24) 


+ 4/x7A(1-^x7A)u3 

The  above  expression  indicates  the  same  order  of  singularity  for  the 
strains  as  for  the  quadrilateral  case. 

In  both  of  the  above  cases  we  have  investigated  the  singularity  at 
Node  1 along  line  1 - 2 of  Figures  1 and  2.  However,  for  the  case  of 
quadrilateral  element,  as  opposed  to  the  case  of  triangular  element, 
the  singularity  at  Node  1 (Figure  1)  along  any  other  ray  emanating  from 
Node  1 is  weaker  than  one-half.  The  singularity  of  order  one-half 
along  any  other  ray  can  be  achieved  by  collapsing  grid  points  1,  4 and 
8 and  placing  grid  points  5 and  7 at  quarter  points  in  cartesian  coor- 
dinates from  grid  point  1. 

Without  loss  of  generality  let  the  cartesian  map  of  such  a col- 
lapsed quadrilateral  element  as  well  as  the  triangular  element  be 
represented  in  Figure  3.  Using  Equations  (7),  (1)  and  (11)  it  can  be 
easily  shown  that 

detlJl  = W(l+c)3sina  (for  the  collapsed  quadrilateral) 

16  2 (25) 

det|Jl  = 2(L2+L3)  sina  (for  the  triangle) 

Both  of  the  above  expressions  vanish  at  grid  point  (1)  along  any  ray. 

The  coordinate  transformations  in  terms  of  cylindrical  coordinates 
r,  e(x  = r cos  e,  y = r sin  e),  are  given  by: 

c = -1  + 2 f cosiera/21  V/2  rl/2 

^ 1 COS  a/2  j 

(For  Quadrilateral  (26) 

= tan(e-a/2)  Element) 

tan  fa/ 2) 


!! 


r1/^2sin(a-8) 

1-2  {2sina  sina/2  cos(e-a/2)  }V2 

r ' sine 

L3  = {2sina  sina/2  cos(o-a/2)}l/2 


(For  Triangular 
Elements) 


When  these  expressions  are  substituted  in  Equations  (2)  and  (13),  it  is 

easily  seen  that  the  displacements  have  the  proper  behavior  to  insure  a 

strain  singularity  of  the  order  one-half  (— ) at  the  crack  tip.  (It 

/r 

should  be  noted  that  the  displacements  for  nodes  1,  4,  and  8 of  the 
collapsed  quadrilateral  must  be  the  same.)  Explicitly, 

u = {l  + fj[-l/4(l-n)(l-f,+n)u2  + l/4(l+n)( C+n-l  )u2 

+ 1/20-0(1-0)^  + l/2(l-n2)^  + l/2(l-0(l+n)u7] 

- orl/2  c_OS1/2(.0-a/2l  I"  0,  COS1/2(e-a/2)  ]/?  / 

2r  COST/2(a/2)  L “ COS^a/Z)  r } \ 


l/4(l-n)u2  - l/4(l+n)u3  + 1/2(1 -n)u5  + l/2(l+n)u7 


- l/4n(l-n)u2  + l/4n(l+n)u3  + 1/2(1  -n  )i^  J 

(For  Quadrilateral 

Element)  (28) 


_ tan(s-a/2) 

n ' tanV/2) 


f'l " 

{2sina  sina/2  cos(e-a/2) }' ^2 
- (u3+4u6)sine J + 0(r) 


(a  +4u^)sin(a-e) 

(For  Triangular 
Elements) 


13 


THE  THREE-DIMENSIONAL  CASE 

The  three-dimensional  twenty-node  isoparametric  quadratic  'Brick' 
element  is  formulated  in  much  the  same  way,  by  mapping  the  geometry 
into  curvilinear  space  U.n.c)  of  a normalized  cube  (-1  * f,  n.cll) 
by  the  quadratic  shape  function  [9], 

20 

x = Z Mc.n.Oxi  , 

i=  i 1 1 

20 

y=l  Ni(L,n,c)yi  , 

i=  i 

20 

z = l Mc.n.cJz*  , 

i=i  (31) 

^ = l/SOHC^O+nniidn^OC^  + nnj  + -2)  l]  n-  q 

+ V4(l-C  2)(l+nn|)(l+CCi)(l-€  •) 

+ l/4(l-n?)(Hffi)(Hui)(l-nJ) 

+ V4(l-c2)(l+«1)(l+nni)(l-cf)  , 

where  N..  is  the  shape  function  at  node  i (i  = 1 to  20)  whose  Cartesian 
and  curvilinear  coordinates  are  (x..,  yi , z.j)  and  , c-)  respec- 

tively'. It  should  be  noted  that  the  shape  function  given  in  Equation 
(31)  is  obtained  by  superposition  of  those  given  in  [9].  The  geometry 
of  the  unit  cube  and  the  numbering  sequence  is  shown  in  Figure  4. 

For  the  isoparametric  formulation,  the  displacements  are  given  by 
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Numbering  sequence  of  an  isoparametric,  quadratic,  'brick' 
element. 


20 

u = E N. (c»n,c)ui  , 

i=  1 

v = z°  Ni(?,n,Ovi  , (32) 

i=  l 

2 0 

w = Z Ni(f,n,rjwi  . 

i=  l 

The  rest  of  the  analysis  follows  in  a similar  fashion  that  given  for  the 
two-dimensional  cases  with  appropriate  augmentation  to  the  three- 
dimensional  quantities. 

The  singularity  element  is  obtained  by  collapsing  one  face,  2376, 
and  placing  the  midside  nodes  9,  13,  11,  15  of  Figure  4 at  quarter 
points.  The  singular  element  is  shown  in  Figure  5,  in  Cartesian  coor- 
dinates. Since  the  elements  are  isoparametric  they  automatically 
satisfy  inter-element  compatibility  and  continuity  in  their  regular  or 
singular  forms.  It  should  be  noted  that  the  displacements  are  not 
singular.  Further,  it  is  easily  shown  that  J N.  = 1.  Hence,  by 
theorems  given  in  [9]  the  elements  satisfy  the  constant  strain  and 
rigid  body  modes.  The  above  conditions  are  necessary  for  the  "patch 
test"  mentioned  in  [2]. 

COMPUTATION  OF  STRESS  INTENSITY  FACTORS 

The  mixed  mode  stress  intensity  factors  can  be  computed  at  quarter 
points  using  the  Westergaard  near  field  displacements  which,  for  plane 
stress,  are  given  in  [10]. 
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Kj  r 1/2 

" “0(2?)  cose/2 
sine/2 


+ sin2e/2 

I +V 


& + «*2«/2 


v = \r_s1/2 


-fc%)'/_  sine/2  |j—  - cos^e/2 
cose/2  [.  |zj+  SI-hVJ 

Solving  for  Kj,  Kjj  we  have: 


Ktt  ’A 

GS' 


Ktt  _ !/2 

+ 

G '21' 


(33) 


K = G(— ) ^2  u cose/2 
I ' r' 


2 v 2 

C - T+J  + cos20/2^  + v sine/2  £t+7  + cos2e/2] 


K11  = G(~) 


l-7+v  - cos2e/2][^7  - cos?o/2] 
2u,  1/2  u sinC/2  - v cos9/2 


(34) 


hfc  ' ccs2e/2] 


It  should  be  noted  that  Equation  (.33)  is  only. a one-term  expansion 
of  the  near  field  crack  solution.  To  obtain  reasonable  accuracy,  the 
quarter  points  have  to  be  very  close  to  the  crack  tip  even  though  the 
elements  surrounding  such  a crack  embody  the  proper  singular  field. 

This  will  increase  the  number  of  elements  required  to  obtain  accurate 
results.  Thus  it  is  advisable  to  use  a higher  order  expansion.  Such 
an  expansion  can  be  obtained  from  airy  stress  functions.  William’s 
stress  function  [11]  is  such  an  expansion  and  using  it,  with  appropriate 
corrections,  the  displacements  under  plane  stress  condition  become. 
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d2n-1  |(n-l/2)cos(n-5/2)e  - 


t 


-1/2 


l 1(-DN  ' 

n=l ,2,3 

JT7  + n - 0 cos(n-l/2)e^j  r" 

(-Dn  d2n  In  cos(n-2)0  - + n + *^ 


cos  ne 


H)"  1a2n_1  (n-l/2)sin(n-5/2)e  + n + 1/2 


(-1 )”a2n  n sin(n_2)0  - + n'j)  sin  ne  ^ 


Comparing  with  the  Westergaard's  solution  we  obtain 


Kj  = - /2ndj 

K„  = ' ^ai 

NASTRAN  IMPLEMENTATION 


(37) 


The  isoparametric  quadratic  quadrilateral,  triangular  and  brick 
elements  have  been  implemented  using  the  NASTRAN  dummy  user  element 
facility  as  outlined  in  Section  6.8.5  of  [12].  This  involved  coding 
element  stiffness  and  stress  data  recovery  subroutines  using  the  analysis 
outlined  above  and  relinking  the  affected  NASTRAN  Links.  Additional 
modifications  were  required  to  some  of  the  Output  File  Processor  (OFP) 
routines.  These  changes  are  detailed  below. 

The  quadrilateral  element  was  implemented  as  a DUM1  element.  KDUM1 , 
the  clc.ent  stiffness  matrix  subroutine,  obtains  material  and  grid  point 
information  from  the  element  connection  and  property  table  (ECPT)  and 
builds  the  matrices  required  to  perform  the  integration  in  Equation  (10). 
The  integration  is  performed  numerically  using  compound  3-point  Gaussian 
quadrature  as  explained  in  [9].  This  results  in  9 evaluations  of  the 
integrand.  Once  the  16  by  16  stiffness  matrix  is  complete,  the  appro- 
priate 2 by  2 submatrices  corresponding  to  the  given  pivot  point  are 
entered  into  the  upper  left  of  the  6 by  6 submatrices  required  by 
SKA1B,  the  stiffness  matrix  insertion  subroutine.  SMA1B  is  called 
8 times  for  each  pivot  point.  The  time  for  element  stiffness  generation 
is  6 seconds  per  element  on  an  IBM  360  Model  44. 


NASTRAN  stress  data  recovery  is  accomplished  in  two  phases.  During 
phase  I,  SDUM11  calculates  [D][B]  from  Equations  (5)  and  (9)  for  each 
grid  point  and  passes  the  resultant  24  by  16  matrix  to  SDUM12  for  final 
stress  calculations.  SDUM11  also  checks  for  singularities  in  the  inverse 
of  the  Jacobian  [Equation  (7)]  and  flags  those  grid  points  which  have  a 
singularity.  Information  passed  to  SDUM12  from  SDUM11  includes  the 
element  id,  grid  point  numbers,  grid  point  singularity  flag,  coordinates 
of  the  eight  grid  points,  crack  orientation  angle  and  the  material  con- 
stants E,  G and  v.  SDUM12,  phase  II  of  the  stress  recovery,  locates 
the  displacements  associated  with  a given  element  and  multiplies  [D][B] 
times  these  displacements  [Equation  (4)  and  (8)]  to  give  the  stress 
components  at  each  of  the  eight  grid  points.  Grid  point  flags  are 
checked  for  singular  grid  points  and,  if  singularities  exist.  Mode  I 
and  Mode  II  stress  intensity  factors  are  calculated  using  Equations  (3b) 
and  (36)  with  a two  term  expansion.  These  stress  intensity  factors  at 
the  quarter  points  are  output  at  the  mid-side  node  of  the  collapsed 
side  while  the  corresponding  corner  node  stress  outputs  are  set  to 
zero.  The  point  ids  of  the  singular  corner  nodes  are  negated  and  the 
mid-side  id  is  set  to  overflow  the  integer  field  specification,  thus 
flagging  the  point  with  asterisks. 

OFP  has  been  modified  to  output  the  eight  sets  of  stress  components 
for  each  element.  These  modifications  were  implemented  by  adding  heading 
formats  to  OFPIA  and  changing  the  appropriate  pointers  and  format 
specifications  in  0FP1BD,  0FP5BD  and  0F1PBD.  Although  the  ADUM  cards 
allow  sufficient  flexibility  to  implement  the  element  stiffness  subroutine. 
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changes  were  required  to  GPTAL  which  describes  the  connection  and 
property  characteristics  of  each  element.  (See  Section  2.5.2.]  in 
[12].)  These  changes  were  required  to  handle  the  expanded  stress 
requirements.  The  number  of  words  SDR2  passes  from  phase  T to  phase 
II  was  changed  from  100  to  430  while  the  count  of  words  SDR2  outputs 
for  real  stresses  was  increased  from  10  to  33. 

The  triangular  element  was  implemented  in  a similar  fashion  to 
the  quadrilateral.  In  this  case  a DUM3  element  was  defined  and 
appropriate  KDUM3,  SDUM31  and  SDUM32  routines  were  written  adhering  to 
the  analysis  defined  in  an  earlier  section.  Numerical  integration  was 
performed  in  KDUM3  using  6 point  Hammer  quadrature  [9]. 

The  implementation  of  the  twenty  node  brick  element  was  not  as 
straightforward  as  that  of  the  quadrilateral  element.  The  brick  was 
implemented  as  a DUM2  element.  DIMENSION  changes  were  required  in 
TA1A  and  TA1B  since  NASTRAN  assumes  a maximum  of  10  grid  points  per 
element.  KDUM2  was  initially  implemented  similarly  to  KDUM1 . However, 
in  this  case,  the  integration  using  a 3-point  Gaussian  quadrature  requires 
27  integrand  evaluations  with  a stiffness  matrix  of  order  60.  The 
size  of  the  KDUM2  subroutine  necessitated  a change  to  the  overlay 
structure  of  LINK3,  placing  KDUM2  in  its  own  overlay  segment.  Also, 
since  NASTRAN  calls  the  stiffness  routines  based  upon  the  pivot  point 
concept,  the  same  brick  element  stiffness  matrix  is  built  twenty  times 
in  an  analysis.  This  technique  resulted  in  a 20  minute  stay  in  SMA1 
for  a one  element  problem  on  an  IBM  360  Model  44. 

KDUM2  was  rewritten  to  build  only  the  requested  submatrices,  t.e. 
those  corresponding  to  the  pivot  point.  Also,  the  matrix  product 


BTDB  (Equation  10)  was  performed  explicitly,  taking  advantage  of  zero 
entries  in  the  B and  D matrices  (Equations  5 and  9).  These  changes 
resulted  in  a reduction  of  computer  time  to  2 minutes  per  element  in 
SMA1 . 

Stress  data  recovery  is  also  nonstandard  for  the  brick  element. 

Due  to  the  size  of  the  arrays  used  in  stress  recovery  for  the  brick, 
phase  I has  a limited  function  of  assembling  arrays  dependent  on 
parameters  not  available  to  phase  II.  The  majority  of  calculations 
for  stress  recovery  are  accomplished  during  phase  II,  saving  storage 
but  increasing  stress  recovery  times.  Quantities  passed  from  phase  I 
to  phase  II  are  the  stress-strain  matrix  [D],  the  element  id,  grid 
point  ids,  grid  point  coordinates  and  the  material  constants  E,  G, 
and  v.  Phase  II  locates  the  displacements,  calculates  the  stress 
components  for  each  grid  point,  checks  for  singular  Jacobians  and 
calculates  stress  intensity  factors  as  required.  Stress  intensity 
factors  are  displayed  in  a manner  similar  to  that  employed  for  the 
quadrilateral  element. 

OFP  has  been  modified  to  output  the  twenty  sets  of  stress  components 
for  each  element.  A DIMENSION  change  was  also  required  in  OFP  to  allow 
twenty  grid  points  per  element.  0FP1A,  0FP1BD,  0FP5BD  and  0F1PBD  were 
updated  to  produce  the  required  heading  and  output  formats.  GPTAL  was 
changed  to  increase  the  number  of  words  to  120  that  SDR2  passes  from 
phase  I to  phase  II.  The  count  of  words  SDR2  outputs  for  real  stresses 
was  increased  to  141. 

The  three  element  implementations  were  checked  independent  of 
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NASTRAN  via  dummy  driver  routines  for  SMA1  and  SDR2.  The  coding  for 
the  element  stiffness  subroutines  was  verified  by  multiplying  the 
stiffness  matrix  for  one  element  times  the  known  displacements  for  a 
uniform  stress  field.  The  resultant  nodal  loading  was  compared  to  that 
found  analytically  from  (for  the  two-dimensional  case) 

{F}e ~_J  1 J1  [B]T  {o}  d£  dn  (38) 

Figure  6 illustrates  equivalent  nodal  forces  for  Oy  = 1 on  the  nor- 
malized square  (-1  <_  x,y  1).  The  stress  recovery  coding  was  checked 
by  passing  the  known  displacements  for  a uniform  stress  field  to  the 
stress  recovery  subroutines  and  observing  the  constant  stress  results. 

The  appendix  describes  the  formats  for  each  of  the  bulk  data  cards. 

A FEW  NUMERICAL  RESULTS 

To  assess  the  accuracy  of  the  method,  three  test  problems  with  a v 
fairly  coarse  grid  (68  elements  and  approximately  239  grid  points)  were 
run.  The  problems  are  the  single  edge  crack,  double  edge  crack  and 
center  crack.  These  three  problems  can  be  done  in  a single  run  as 
subcases  with  different  single  point  constraints  as  shown  in  Figure 
7a,  b,  and  c.  The  solutions  of  the  above  problems,  by  various  methods, 
have  been  well  documented  [13].  Using  quadrilaterial  elements,  it 
was  found  that  the  finite  element  solutions  were  accurate  to  within 
2-3%.  Graphically,  this  is  illustrated  in  Figure  8 for  the  double 
edge  crack  by  using  the  Westergaard  near  field  solution  [10]  for  Oy 

ASTM  has  stringent  requirements  for  the  size  of  specimens  for 
fracture  toughness  testing.  However,  in  many  applications  of 
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• CENTER  CRACK 


thick-walled  cylinders  these  requirements  are  not  easily  met.  The 
C-shaped  specimen,  which  is  easily  obtained  from  thick-walled  cylinders, 
was  suggested  [14]  and  is  now  accepted  as  a standard  test  for  such 
cylindrical  material.  The  stress  intensity  factors  for  such  a section, 

shown  in  Figure 9,  were  computed  for  different  crack  lengths  and  the 

finite  element  results,  experimental  results,  and  the  collocation 
results  of  [15]  are  shown  in  Figure  9 . It  is  seen  that  remarkable 
agreement  is  obtained  with  just  48  elements  and  171  grid  points.  The 
results  have  also  been  compared  with  those  in  [16]  and  similar  corre- 
spondence was  observed. 

THE  STABILITY  OF  QUADRILATERAL  AND  TRIANGULAR  ELEMENTS 

During  the  process  of  implementation  of  these  elements  into 
NASTRAN  it  was  found  that  a slight  perturbation  in  placing  the  mid- 
side node  opposite  to  the  crack  tip  for  a collapsed  quadrilateral 
element  led  to  a substantial  error  in  the  computation  of  stress  inten- 
sity factors.  (See  the  Conclusion  in  [18].)  No  such  errors  could  be 
detected  for  triangular  elements.  In  this  section  it  is  demonstrated 
that  under  perturbation  the  triangular  element  preserves  the  singularity 
required  for  crack  problems  while  the  collapsed  quadrilateral  element 
does  not  [19]. 

TRIANGULAR  ELEMENT 

Let  the  grid  point  5 be  perturbed  as  shown  in  Figure  10a. 

(Compare  with  Figure  3.)  Denoting  the  perturbed  quantities  with  an 
asterisk  we  have 
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Figure  9.  Comparison  of  NASTRAN  results  with  those  obtained 
experimentally  and  by  collocation. 
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We  have 


x5  = cos  ~ + a 


y5  = sirv|  cos-|  + b 


For  Grid  Point  5 


Det | J* | = Det|J|  - (8l_2sina)a  - (8L2cos<i  - 4L0)b 


As  seen  the  perturbed  Jacobian  vanishes  at  grid  point  1 (Lj  = 1, 
L2  = L3  = 0)  for  all  values  of  perturbation  parameters  a and  b. 

order  of  singularity,  however,  will  only  become  apparent  when  we 
at  the  coordinate  transformation 


x -.E  = x + aNc 

1=  l 1 1 5 


= (L2  + L3)[L2  + L3cosa]  + 4aL2L3 
y*  j Niyi  = y + bN5 

= (L2  + L 3 ) L3  s i net  + 4bl.2L3 


Solving  for  L,  we  have 


L2  _ (By  + cx*)  * /(By^  + cx*)7~-4Ay*7 

3 ~ 2A 

where 

2 2 
A = -4a  sin  a - 4b  sina(l-cosa)  - 16ab  sina  + 16b  cosa 

B = 2sina  - (4b  + sina)(l  + 4a  + cosa)  (4 

C = (4b  + sina)2 

Using  the  polar  coordinates  (x*  = r coso,  y*  = r sino)  we  have  from 
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= r1/2  (Bsing  + 


sing  + Ccose)2  - 4As 
2A 


in2e  j 


with  a similar  expression  for  L2.  Comparing  Equation  (45)  with  the 
unperturbed  case  [Equation  (27)]  we  see  that  L2,  L3  and  consequently 
the  displacements  have  the  same  behavior  as  the  unperturbed  case. 
Hence  for  the  case  of  the  triangular  element,  the  singularity  is  pre- 
served under  perturbation. 

QUADRILATERAL  CASE 

Let  grid  point  6 be  perturbed  as  shewn  in  Figure  10b.  For 
simplicity  let  a = b = e.  (It  is  sufficient  to  show  the  lack  of 
singularity  on  one  ray.)  Then 

x*  = l(l+c)2  ^(1+cosa)  - n(l-C0Sa)J  + |_(1  -n2 ) (1  +0  ( 

y*  = lsina(l+c)2(l+n)  + f(l-n2)(l  + 0 ( 

Along  the  line  x*  = y*  we  have  from  (46)  and  (47), 

n = j + cos-a--~  -slna  = constant  = c (SAY) 

1 - cosa  + sina 

Substituting  n into  Equation  (47)  and  solving  we  have  (in  terms  of 
polar  coordinates,  x*  = y*  = r sine,  e = u/4) 


n + 0 = 1 ' r Sine  (48) 

Since  (1+4)  is  the  common  factor  in  the  displacement  (28)  it  is  seen 
that  there  is  at  least  one  ray  in  the  collapsed  quadrilateral  case  where 
the  singularity  required  for  the  crack  problems  disappears.  Hence  such 
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e.'iurbations  of  the  collapsed  quadrilateral  elements,  which  may  either 
be  caused  by  automatic  mesh  generation  or  simply  by  transcription  errors 
may  lead  to  unstable  answers. 

NUMERICAL  RESULTS  FOR  STABILITY 

To  illustrate  the  effect  of  perturbing  the  mid-side  node  let  us 
consider  a crack  of  length  2a  in  an  infinite  media  under  a uniform 
tension  of  P applied  normal  to  the  crack.  The  solution  of  the  problem 
is  well  known.  Using  Muskhel ishvili ‘s  method  [17],  we  have  for  the  plane 
stress  condition  displacements  given  by 

u + iv  = ^ ffej  - 2 - ^TT  J (49) 


where 


v ( Z ) = P* 


* f-  1 


l^(z)  = 


2 z/a 


(£  - ') 


[f:r7- 


1/2 ' 


For  a finite  element  grid  of  12  elements  as  shown  in  Figure  11,  the 

displacements  computed  from  (49)  were  prescribed  at  corner  grid  points. 

The  displacements  were  computed  at  quarter  points  using  collapsed 

quadrilateral  as  well  as  triangular  elements.  The  results  obtained 

by  both  methods  were  in  very  close  agreement  with  the  theoretical 

results  of  (49).  However  when  the  opposite  node  of  the  1st  element 

was  perturbed  the  results  from  the  collapsed  quadrilateral  elements 

were  quite  unstable  compared  to  those  obtained  from  triangular 

elements.  This  is  illustrated  in  Figure  (12)  with  0%,  20%,  and 

40%  perturbation  of  a single  opposite  node  of  the  1st  element.  It 

was  seen  that  the  error  is  of  the  order  — for  the  quadrilateral 

r 

case.  There  was  no  significant  change  in  results  for  the  case  of 
triangular  elements  as  seen  in  Figure  13. 

CONCLUSIONS 

Quadratic  quadrilateral  and  triangular  isoparametric  elements  can 
be  easily  used  to  construct  an  elastic  singular  element  for  analyzing 
crack  problems.  Collapsed  quadrilateral  elements  are  less  stable  in 
preserving  the  elastic  singularity,  under  perturbation  of  the  opposite 
mid-side  nodes,  than  the  triangular  elements.  These  elements  have  been 
implemented  in  NASTRAN.  Similar  elements  can  be  also  constructed  from 
cubic  elements.  The  plastic  singularity  can  also  be  achieved  by  the 
use  of  "sliding  node"  concepts.  These  questions  will  be  investigated 
in  the  near  future. 
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Finite  element  grid  to  test  the  effect  of  perturbation  of  a node,  with  displacements 
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BULK  DATA  DECK 


Input  Data  Card  ADUM1 

Description:  Defines  attributes  for  an  isoparametric,  quadratic  quadrilateral 
implemented  as  a DUM1  dummy  user  element. 

Format  and  E xample: 


ADUM1 

NG 

NC 

NP  ND 

ADUM1 

8 

0 

1 3 

Field  Contents 


NG  Number  of  grid  points  connected  by  DUM1  (Integer  = 8) 

NC  Number  of  additional  entries  on  CDUM1  (Integer  = 0) 

NP  Number  of  additional  entries  on  PDUM1  (Integer  = 1) 

ND  Number  of  displacement  components  at  each  grid  point 

used  in  generation  of  differential  stiffness 
(Integer  = 3) 

Remarks : 

1.  All  fields  must  correspond  to  the  example. 

2.  One  ADUM1  must  be  present  if  any  CDUMl's  are  present. 


BULK  DATA  DECK 


Input  Data  Card  ADUM3 

Description:  Defines  attributes  for  an  isoparametric,  quadratic  triangle 
implemented  as  a DUM3  dummy  user  element. 

Formal  and  Example^ 


ADUM3 

NG 

NC 

NP 

ND 

ADUM3 

6 

0 

1 

3 

Field  Contents 


NG  Number  of  grid  points  connected  by  DUM3  (Integer  = 6) 

NC  Number  of  additional  entries  on  CDUM3  (Integer  = 0) 

NP  Number  of  additional  entries  on  PDUM3  (Integer  = 1) 

ND  Number  of  displacement  components  at  each  grid  point 

used  in  generation  of  differential  stiffness 
(Integer  = 3) 

Remarks: 

1.  All  fields  must  correspond  to  the  example. 

2.  One  ADUM3  must  be  present  if  any  CDUM3's  are  present. 


■ 
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BULK  DATA  DECK 
Input  Data  Card  ADUM2 

Description:  Defines  attributes  for  an  isoparametric,  quadratic  brick 
implemented  as  a DUM2  dummy  user  element. 

Format  and  Example: 


ADUM2 

NG 

NC 

NP 

ND 

ADUM2 

20 

0 

1 

3 

Field  Contents 

NG  Number  of  grid  points  connected  by  DUM2  (Integer  = 20) 

NC  Number  of  additional  entries  on  CDUM2  (Integer  = 0) 

NP  Number  of  additional  entries  on  PDUM2  (Integer  = 1) 

ND  Number  of  displacement  components  at  each  grid  point 

used  in  generation  of  differential  stiffness 
(Integer  = 3) 

Remarks: 

1.  All  fields  must  correspond  to  the  example. 

2.  Cne  ADUM2  must  be  present  if  any  CDUM2’s  are  present. 
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Remarks: 


1.  The  order  of  the  grid  points  is:  Cl,  C2,  C3,  C4  counterclockwise; 
Ml  is  between  Cl  and  C2,  followed  by  M2,  M3,  M4. 

2.  The  quadrilateral  must  lie  wholly  in  the  X-Y  plane. 

3.  Stresses  are  given  in  the  global  system.  For  collapsed  elements, 
mode  I and  II  stress  intensities  are  output. 

4.  If  one  side  is  collapsed  to  form  a crack  element,  these  rules  must 
be  followed: 

a.  Sides  adjacent  to  the  collapsed  side  must  be  straight  and  have 
their  'mid-side1  nodes  located  exactly  1/4  of  the  length  of  the 
side  from  the  collapsed  side. 


44 


mmm 


b.  The  side  opposite  the  collapsed  side  must  be  straight  and  have 
its  'mid-side'  node  located  exactly  1/2  the  distance  between  its 
corresponding  corner  nodes. 

5.  If  degrees  of  freedom  on  the  collapsed  side  are  not  constrained  via 
single-point  constraints,  they  should  be  constrained  via  multi- 
point constraints  to  maintain  a physical  correspondence  (i.e., 

Ui  = U4  = U8,  Vi  = V4  = V8). 
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BULK  DATA  DECK 


Input  Data  Card  CDUM2 

Description:  Defines  an  isoparametric,  quadratic  bricP  element  (DUM2), 
Format  and  Example: 


Field 

EID 

PID 

Cl , . . ,C8 
Ml,.., Ml  2 


CDUM2 

EID 

PID 

Cl  C2 

C3 

C4 

C5 

C6 

+abc 

CDUM2 

31 

7 

1 2 

3 

4 

5 

6 

ABC 

+abc 

C7 

C8 

Ml 

M2 

M3 

M4 

M5 

M6 

+def 

+BC 

7 

8 

9 

10 

11 

12 

13 

14 

DEF 

+def 

M7 

M8 

M9 

Ml  0 

Mil 

Ml  2 

+EF 

15 

16 

17 

18 

19 

20 

Contents 

Element  identification  number  (Integer  > 0) 
Identification  number  of  a PDUM2  property  card 
(Integer  > 0) 

Grid  point  identification  of  corner  points  (Integers 
> 0,  Cl  t...f  C8) 

Grid  point  identification  numbers  of  'mid-side' 
points  (Integers  > 0,  Ml  ?...  ? M12) 


Remarks: 


1.  The  order  of  the  grid  points  is: 

Cl,  C2,  C3,  C4  about  one  face 

C5,  C6,  C7,  C8  about  the  opposite  face 

Ml  is  between  Cl  and  C2,  followed  by  M2,  M3,  M4 

M5  is  between  C5  and  C6,  followed  by  M6,  M7,  M8 

M9  is  between  Cl  and  C5,  followed  by  M10,  Mil,  Ml 2. 

2.  Stresses  are  given  in  the  global  system.  For  collapsed  bricks, 
mode  I,  II  and  III  stress  intensities  are  given. 

3.  For  collapsed  elements,  the  crack  must  lie  wholly  within  some 
X-Y  plane. 

4.  If  one  face  is  collapsed  to  form  a crack  element,  these  rules  must 
be  followed: 

a.  Sides  adjacent  to  the  collapsed  face  must  be  straight  and  have 
their  'mid-side'  nodes  exactly  1/4  of  the  length  of  the  side  from 
the  collapsed  face. 

b.  The  line  formed  by  the  collapse  of  a face  need  not  be  straight. 

5.  If  degrees  of  freedom  on  the  collapsed  side  are  not  constrained  via 
single-point  constraints,  they  should  be  constrained  via  multi- 
point constraints  to  maintain  physical  correspondence  (i.e., 

Ui  = (Jit  = U5,  V)  = V17  = V 5,  Wi  = Wn  = W5,  U9  = U] 3 , V 9 = V 1 3 , 

etc. ). 
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BULK  DATA  DECK 
Input  Data  Card  CDUM3 

Description:  Defines  an  isoparametric,  quadratic,  triangular  element 
(DUM3) 

Format  and  Example: 


CDUM3 

EID 

PID 

Cl  C2 

C3 

Ml 

M2 

M3 

CDUM3 

25 

14 

1 2 

3 

4 

5 

6 

Field 


EID 

PID 

Cl ,C2  ,C3 
Ml ,M2,M3 


Remar  ks_: 

1.  The  order  of  the  grid  points  is:  Cl,  C2,  C3  counterclockwise;  Ml  is 
between  Cl  and  C2,  followed  by  M2,  M3. 

2.  The  triangle  must  lie  wholly  in  the  X-Y  plane. 

3.  Stresses  are  given  in  the  global  system.  For  singular  elements, 
mode  I and  II  stress  intensities  are  output. 

4.  To  form  a crack  element,  these  rules  must  be  followed: 

a.  Sides  adjacent  to  the  crack  tip  must  be  straight  and  have  their 
Vnid-side'  nodes  located  exactly  1/4  of  the  length  of  the  side 
from  the  crack  tip. 

b.  The  side  opposite  the  crack  tip  must  be  straight  and  have  its 
'mid-side'  node  located  exactly  1/2  the  distance  between  its 
corresponding  corner  nodes. 


Contents 

Element  identification  number  (Integer  > 0) 
Identification  number  of  a PDUM1  property  card 
(Integer  > 0) 

Grid  point  identification  numbers  of  corner  points 
(Integers  > 0,  Cl  / C2  / C3) 

Grid  point  identification  numbers  of  'mid-side' 
points  (Integers  > 0,  Ml  / M2  / M3) 
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BULK  DATA  DECK 


Input  Data  Card  PDUM1_ 

Description:  Defines  the  properties  of  an  isoparametric,  quadratic, 
quadrilateral  element.  Referenced  by  a CDUM1  card. 

Format  and  Example: 


PDUM1 

PID 

MID 

BETA 

PDUM1 

23 

5 

45.0 

Field  Contents 

PI D Property  identification  number 

MID  Material  identification  number  of  a MAT1  card 

BETA  Orientation  angle  in  degrees  of  crack  with  respect 

to  global  system  (Blank  or  real).  BETA  is  positive 
for  the  top  elements  of  the  crack  and  negative  for 
the  bottom  elements  of  the  crack. 


Remarks: 

1.  If  RfTA  is  blank,  0.0.  is  assumed. 

2.  SETA  is  only  used  for  elements  with  a collapsed  side.  If  specified 
for  other  elements,  it  is  ignored. 
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BULK  DATA  DECK 


Input  Data  Card  PDUM2 

Description:  Defines  the  properties  of  an  isoparametric,  quadratic, 
brick  element.  Referenced  by  a CDUM2  card. 

Format  and  Example: 


PDUM2 

PID 

MID 

PDUM2 

75 

1 

Field 

Contents 

PID 

MID 

Property  identification 
Material  identification 

number 

number  of  a MAT1  card 

BULK  DATA  DECK 


Input  Da^a  Card 


PDUM3 


Description:  Defines  the  properties  of  an  isoparametric,  quadratic, 
triangular  element.  Referenced  by  a CDUM3  card. 


Format  and  Example: 


PDUM3 

PID 

MID 

BETA 

PDUM3 

23 

5 

45.0 

Field  Contents 

PID  Property  identification  number 

MID  Material  identification  number  of  a MAT1  card 

BETA  Orientation  angle  in  degrees  of  crack  with  respect  to 

global  system  (blank  or  real).  BETA  is  positive 
for  the  top  elements  of  the  crack  and  negative  for 
the  bottom  elements  of  the  crack. 


Remarks : 

1.  If  BETA  is  blank,  0.0.  is  assumed. 

2.  BETA  is  only  used  for  elements  with  a singularity.  If  specified 
for  other  elements,  it  is  ignored. 
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